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ABSTRACT. Incompressible, homogeneous magnetohydrodynamic (MHD) turbulence 
consists of fluctuating vorticity and magnetic fields, which are represented in terms of 
their Fourier coefficients. Here, a set of five Fourier spectral transform method numerical 
simulations of two-dimensional (2-D) MHD turbulence on a 5 1 2 2 grid is described. Each 
simulation is a numerically realized dynamical system consisting of Fourier modes 
associated with wave vectors k, with integer components, such that k = |k| < & max . The 
simulation set consists of one ideal (non-dissipative) case and four real (dissipative) 
cases. All five runs had equivalent initial conditions. The dimensions of the dynamical 
systems associated with these cases are the numbers of independent real and imaginary 
parts of the Fourier modes. The ideal simulation has a dimension of 366104, while each 
real simulation has a dimension of 411712. The real runs vary in magnetic Prandtl 
number Pm, with {0.1, 0.25, 1, 4}. In the results presented here, all runs have been 
taken to a simulation time of t = 25. Although ideal and real Fourier spectra are quite 
different at high k, they are similar at low values of k. Their low k behavior indicates the 
existence of broken symmetry and coherent structure in real MHD turbulence, similar to 
what exists in ideal MHD turbulence. The value of Pm strongly affects the ratio of kinetic 
to magnetic energy and energy dissipation (which is mostly ohmic). The relevance of 
these results to 3-D Navier-Stokes and MHD turbulence is discussed. 
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1. Introduction. Sufficiently energetic continuum flow, whether it consists of liquid, 
gas, or plasma, tends to develop instabilities that lead to the large-scale chaotic motion 
called turbulence. Although fluid compression may be present initially, density 
fluctuations can radiate away as sound waves or damp out due to bulk viscosity, leaving a 
large-scale flow that is essentially incompressible. This incompressible limit is 
characterized by its fluid vorticity, and additionally, in the case of plasma, by its 
magnetic field. Understanding magnetohydrodynamic (MHD) turbulence requires solving 
the equations describing the evolution of the vorticity and magnetic induction. However, 
due to the inherent non-linearity of these equations, there are no general analytical 
solutions (as linear equations have) and they must be solved numerically. 

A fundamental focus for research is in an area termed homogeneous turbulence [1], 
where the region of fluid under consideration is assumed to be remote from any physical 
boundaries. In this case, the appropriate numerical boundaries are periodic and Fourier 

series are used as the basis for analytical approaches and computational techniques. 
Fourier transformation turns the few basic partial differential equations describing the 
flow in physical space (x-space) into many ordinary differential equations (ODEs) in the 
Fourier transform space (A>space). To enable numerical solution, these ODEs must be 
finite in number, requiring that finite Fourier series be used. Typically, the truncation is 
done by retaining only those Fourier modes whose wave vectors k (with integer 
components) have magnitudes less than or equal to some maximum value, k - |k| < k max 
(for dissipative turbulence k max = ko, the dissipation wave number). The number of 
independent real and imaginary parts of the Fourier coefficients within the ball k < k max 
determines the dimension of the dynamical system embodied by these &-space ODEs. 

Since each discrete Fourier mode in £-space represents a continuous function in x-space, 
we have a discrete dynamical system (the Fourier modes) representing, to a certain level 
of resolution, a continuous physical system. Although the Fourier modes are themselves 
continuous analytical functions of time, the actual values they can take are discrete, being 
drawn from the possible digital states of a computer word (of 64 bits in the numerical 
simulations to be described herein). Thus we have, in computational physics, the 
unavoidable result that continuum systems are always represented by discrete dynamical 
systems, consisting, as they do, of algorithms, their numerical implementation, and the 
digital hardware on which they are run. The viability of these simulations is ultimately 
determined through comparison with theoretical prediction and experimental 
measurement. 

Here, we will concentrate our efforts on examining certain aspects of two-dimensional, 
homogeneous, magnetohydrodynamic turbulence (hereafter simply called 2-D MHD 
turbulence). In particular, we will consider both ideal (non-dissipative) and real 
(dissipative) 2-D MHD turbulence. In the ideal case, there is a statistical theory with 
which to compare numerical results [7, 10, 16, 17, 30-37], and in the real cases we 
compare the distribution of energy with inertial decay laws [14, 24] and with the 
distribution found in the ideal case. However, we do not add energy through spectral 
forcing to the real case, so that the dissipative simulations are autonomous dynamical 
systems [38] in which energy decays with time. 
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The one ideal and four real cases presented here all start from equivalent initial 
conditions, so that a controlled comparison can be made. The real cases differ in their 
values of kinetic viscosity v and magnetic diffusivity r| (so that the runs are identified by 
their magnetic Prandtl number. Pm = v/r|). These five runs took about three months of 
cpu-time on a DEC Alpha, allowing each simulation (with 64-bit words) to be taken up to 
a simulation time of t = 25. (In the future, we plan to extend these times by an order of 
magnitude to study long-term states of decaying 2-D MHD turbulence with varying Pm', 
the runs to be presented contain significant new results that motivate this extension). 

The growth and quasi-stationarity of low -k modal values in decaying 2-D MHD 
turbulence, similar to the behavior seen in the ideal case, is a primary result. In particular, 
in the real runs those Fourier modes with low values of k exhibited manifestations of the 
dynamically broken symmetry and associated coherent structure that has been found in 
ideal turbulence [32, 33, 35]. Thus, although the ideal and real Fourier spectra are quite 
different at high k, their similarity at \ow-k, indicates an area of overlap between ideal 
and real MHD turbulence. 


These results are also relevant to 3-D homogeneous Navier-Stokes and MHD turbulence. 
The reason is that homogeneous 2-D and 3-D MHD turbulence, as well as 3-D Navier- 
Stokes, each have, in addition to energy, an ideal integral invariant called a helicity [2, 6, 
23, 39]. It is the presence of helicity that dynamically breaks the symmetry of the basic 
equations[32, 33, 35]. (The ideal invariants for homogeneous 2-D Navier-Stokes 
turbulence are energy and enstrophy, neither of which is helical.) In this regard, 2-D 
MHD turbulence serves as representative of 3-D turbulence, one that allows a 
substantially greater k mm than possible in 3-D simulations. 

In outline, the following sections contain a brief, but thorough discussion of the basic 
MHD equations, their 2-D form, and their 2-D integral invariants. Fourier transforms are 
then discussed, along with some aspects of 2-D MHD turbulence as a dynamical system. 
Next, the statistical mechanics of ideal 2-D MHD turbulence, the spectral laws for real 
turbulence, the numerical method and results, and finally, concluding remarks are 
presented. 

2. Basic equations. The first of the basic equations of MHD turbulence is the Navier- 
Stokes equation with the electromagnetic force term jxB [20, 28] (where j is the electric 
current and B is the magnetic induction): 


on „ 

1- u • Vu 

dt 


= -V/? + pV u + jxB 


( 1 ) 


Here, p is the mass density,/? is the pressure and p is the dynamic viscosity. The dynamic 
viscosity is related to the kinematic viscosity v through p = pv. The assumption of 
incompressibility means that p = p 0 = constant or equivalently that Vu = 0. 
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The magnetic evolution equation originates in Maxwell’s equations [12], which in SI 
(Systeme International) units (i.e., the MKSA system), are 

a) — = -VxE b) V • B = 0 (2) 

dt 

a) ^ = VxH-j b) V-E = p e (3) 

dt 

Above we have the electric field E, the electric displacement D, the magnetic field H and 
the magnetic induction B. We assume that the fluid is neither dielectric nor permeable, so 
that the electric variables E and D are related through the free space electric permittivity 
e 0 , and the magnetic variables H and B are related through the free space magnetic 
permeability p 0 (both e 0 and p 0 are constants), while a simple form of Ohm’s law 
connects j and E (here, the electrical conductivity a is assumed constant): 

a) D = e 0 E b) B = p 0 H c) j = aE. (4) 

In MHD, incompressibility implies that velocities are small with respect to the speed of 
sound and therefore also to the speed of light. Thus, electromagnetic radiation is ignored, 
so that on the left side of (3a) the displacement current 0D/8/ is set to zero. (Since the 
addition of the displacement current was Maxwell’s unifying, in dropping it, (2) and (3) 
become the ‘no yet Maxwell equations’ [12]) 

Neglecting the displacement current in (3a) and using (4b) gives the relation between 
electric current and magnetic induction: 

Moj= Vx B. (5) 

Obviously, V j = 0, so that both B and j have zero-divergence. 

Now, combining Faraday’s law of magnetic induction (2a), Ohm's law (4c), and the 
current-magnetic induction relation (5), yields the second basic equation of MHD, the 
magnetic induction equation : 

^ = Vx(uxB) + r|V 2 B. (6) 

1 *2 
The magnetic diffusivity r| = (|T 0 o)' is also a constant ( r\ and v have the same units: m Is 

in SI units). In MHD, equation (6) determines B and 'B determines j’ through equation 

(5), rather than the usual ‘j determines B' of Maxwellian electrodynamics. 

Equations (1) and (6) are the basic equations of MHD. Since incompressibility is 
assumed, the system of equations can be non-dimensionalized, so that the constant p is 
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eliminated from (1). The pressure p is removed from consideration by taking the curl of 
equation (1) (with p = 1) to yield an equation for the vorticity 0) = Vxu of the magneto- 
fluid: 


— = Vx(uxco + jxB) + vV 2 o> (7) 

dt V ’ 

The pressure, if required, can be found by taking the divergence of (1). 

Equations (6) and (7) are the basic equations of MHD turbulence. In (7), the inverse of v 
is the Reynolds number, while in (6) the inverse of r| is the magnetic Reynolds number. 
Note that all of the fields u, co, B and j appearing in (6) and (7) have zero divergence. In 
the case of B, this leads to the well known magnetic vector potential A, i.e., B = VxA. 

3. 2-D MHD. In the presence of a constant, externally imposed uniform magnetic 
induction B 0 , the total induction B can be split into B 0 and a fluctuating part b: 

a) B = B 0 + b b) Vb = 0. (8) 

In 3-D cartesian coordinates x, y, z, the corresponding orthonormal basis vectors are e x , 
e y , e z , respectively. Let us choose B 0 = B 0 e z . In the limit B 0 — > °°, with fixed limits on 
fluctuating magnetic and kinetic energies, equations (6) and (7) give 


a) 


T 

„ a 

|W| 

= B n — 



dz 



b) ±(a = ] = f(x,y,z±V A t). 


( 9 ) 


The Alfven velocity is V A = B 0 in these dimensionless equations. If the solutions (9b) 
represent initially localized disturbances in an infinite medium, they propagate away, 
leaving a 2-D flow field. In the case of a bounded flow, such as homogeneous MHD 
turbulence, there is nowhere to go. However, in physical experiments [40] and in 
numerical simulations [30, 31, 25], MHD turbulence has been seen to lose its z 
dependence for increasing B 0 . Thus, in the limit of large B 0 , 3-D MHD turbulence 
becomes effectively 2-D, and the use of 2-D simulations to represent physical phenomena 
becomes viable. 


Also, a 2-D approach is advantageous numerically, since a denser grid can be used than is 
the case in 3-D. For example, if we use a 512x512 grid in 2-D, we have 512 points per 
dimension, or 2 18 points overall. In 3-D, the same number of points would be distributed 
evenly into a 64x64x64 grid, but with only Vs as many points along each dimension. 
Hence, restricting MHD to 2-D, in addition to being appropriate in the presence of a very 
large mean magnetic field, also allows for a greatly increased spatial resolution. For these 
reasons, we now develop and then numerically solve the 2-D MHD equations. 

Since V b = V j = V u = Vco = 0, we can define for 2-D MHD the following: 
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a) y/= ip(x,y,t)e z 

b) a = a(x,y,t)e z 

(10) 

a) u = Vxip- - e z xV ip 

b) b = Vxa = - e z xV«. 

(11) 


The quantity y/(x,y,t) is the stream function and a(x,y,t ) is the magnetic scalar potential. 
In terms of these, the vorticity and current are 


a) to = Vxu = coe z 

b) j = Vxb =je z 

(12) 

a) 0) = -VV 

b) j = -V 2 o. 

(13) 


If we now put the expressions (10) through (13) into (7) and (6), we get 

— +uVco = b-Vy + vV 2 0) (14) 

dt 


— + u-Va = riV 2 o. (15) 

dt 


In these equations, the spatial derivatives are only taken with respect to x and y, e.g., 


u Vo = 


da dip 
dx dy 


da dip 
dy dx 


J(a,ip) 


(16) 


V 


a = 


d 2 a d 2 a 

al 2 ’ + a7"' 


(17) 


The basic equations of 2-D MHD turbulence are (15) and (16). Since we assume 
homogeneity, the various scalar fields a, \p, j and co are expanded in terms of truncated 
Fourier series, and Fourier transformation takes the two partial differential equations (14) 
and (15) into a great number of ordinary differential equations. We will move on to this 
topic in the next section. 

First, however, we will demonstrate the existence of integral invariants in a periodic 
domain in x-space. It is straightforward to show, using (11a), that the Navier-Stokes 
equation (1) (with p = 1) can be written as 

^ = -(oe : xV^-e,xv(^ + 4 u 2 )+ je z xVa-vVty . (18) 


Taking the gradient of (16) gives a similar equation for Vo: 
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dVa 

dt 


= J(Vi//,a) + J(i//,Va)-r\VJ 


(19) 


At this point, we take the dot product of (18) with Vt/t, and of (19) with Va, and adding 
the results, to get the energy flux equation: 


0 £ 

a) — - = -V • (j? w u + (u • Vcr)Va]- vVi// ■ Vco-r|Va • VJ 

dt 


( 20 ) 


a) £— 4-(|u| 2 + |b| 2 ) 


b) P u = P + \ u2 


( 21 ) 


Next, a times (15) plus oo times ( 1 6) gives 


d— + V -(acou - ajb) = vaV 2 co-r|CD/ . 

dt 


( 22 ) 


Lastly, a times (16) yields 

a « 2 


dt 


+ V - (a u) = -2r\aj 


(23) 


Now, we integrate (20), (22) and (23) over a periodic square of side length 2n, to find the 
volume averages of the quantities e, aw and a 2 , respectively. The volume average of any 
quantity q will be defined as 

- 2n 2n 

[^r] s (2n) 2 \dy\dx q(x, y) . (24) 

o o 

All of the divergence terms will vanish because of the periodic boundary condition and 
the terms containing v and r) can be integrated by parts. The results are 


dE 

— = -2(vQ + ri J) 

dt 


(25) 

a) E = [ e ] b) Q. = 

l[ro 2 ] c)J = ![/]. 

(26) 

C = (v + Y\)Hj 
at 


(27) 

a) H c =-j[aco] 

b) Hj =l[ya)] 

(28) 
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( 29 ) 


a) A = y[a 2 ] 


b) E m = ±[\b\ 2 ]. 


(30) 


For completeness, we also define 


a) E k ={[|u| 2 ] 


b )E = E k +E m . 


( 31 ) 


In the above equations, E is the total mean energy, while Ek and Em are the mean kinetic 
and magnetic energies, respectively. The mean squared vorticity is Q. (called the 
enstrophy), while J and A are the mean squared current and mean squared magnetic 
potential , respectively. The integral He is called the cross helicity, while Hj is also a 
helicity (/. e . , a pseudoscalar integral), though it is not an ideal invariant. 

If we set the viscosity and diffusivity equal to zero, v = r| = 0, we have ideal MHD. In 
this case, equations (25), (27) and (29) tell us that E, He and A are integral invariants. In 
fact, these are the only ideal invariants in 2-D MHD (and they remain so when the fields 
a and w are represented by truncated Fourier series, as will be demonstrated presently). If 
v ^ 0 and r| ^ 0, we have real MHD, wherein E and A decay monotonically, while He 
may increase or decrease. In what follows, we present results from numerical simulations 
of both ideal and real MHD turbulence. First, however, we discuss the transformation to 
£-space and the explicit dynamical systems that ensue. 

4. Fourier transform methods. In x-space, 2-D MHD turbulence is represented by the 
value of its fields a(x ) and oo(x) (and the fields analytically related to these) within an 
area of side length 2n, where the position vector is x = xe x + ye y , with 0 <x,y < 2n. Since 
o(x) and co(x) are assumed to be periodic in this area they will be represented by discrete 
Fourier series, e.g 


Note that time t is implicit in the arguments of all variables such as a(x) and d(k), being 
omitted only for brevity. Also, N = 2 M for use in fast-Fourier transforms (FFTs) [4]. Here, 
we use 2-D wave vectors k, with integer components k x and k y : 


a) co(x) = 1 Ll co(k)e /kx b) w(k) S(x)<T' k x . (32) 



a) k = k x e x + k y e y 



(33) 


The x-space position vectors x is transform (32b) are also discrete: 


a) x = -2a- (ne x + me y ) 


b) 0 <n,m<N. 


(34) 
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The x-space values (o(x) are real, so the Fourier coefficients satisfy the ‘reality condition’ 
* * 

(o(k) = co (-k), where denotes complex conjugation[henceforth we drop the ~ over c/(k), 
co(k), etc.]. Thus, (32), (33) and (34) define an invertible transformation between N 
values ct)(x) and the N values Re co(k), lm co(k) [and similarly for a(k)]. 

In a numerical implementation, we wish to maintain isotropy in &-space and to avoid 
aliasing effects in the evaluation of the non-linear terms in (15) and (16). Thus, we 
require that |k| < k max < N/2, which imposes what is termed isotropic truncation in k- 
space for an N-by-N transformation. For a real (dissipative) non-linear simulations, Ar max 
is chosen as (A/2) 2 - !4; in these pseudospectral simulations [3, 5, 9], dissipation renders 
aliasing effects negligible. In the case of an ideal simulation, a de-aliasing procedure is 
required to evaluate non-linear terms; we use the method due to Patterson and Orszag 
[27] , in which * ma x = 2N 2 /9, to create a Fourier spectral transform method [3, 5, 9]. 

Fourier transformation of the non-linear evolution equations (15) and (16) produces the 
following two equations in £-space: 

^ k ~ = i I e ; pxq[^(p)(u(q) + y(p)o(q)]-v/: 2 co(k) (35) 

dt p+q=k 


^l = i S e z • pxq a(p)(c/(q) -r|& 2 c/(k) . (36) 

dt p+q=k 

If we multiply (35) by yf (k) and (36) by j*{ k), and sum over k (all sums are such that |k|, 
|p|, |ql <^max) we have 

dE 

— = 6(y/, if/, co) + 6(\\i, j, a) + 6(a, j, y) - 2 (vQ + r(J) (37) 

dt 

6(a,J,<f/) = -i Z e. -pxqa(k)/(p)y/(q) • (38) 

p+q+k=0 

In the triple summation in (38), k. p and q are dummy summation variables and can be 
interchanged to yield 

9(a, ;>) = -0(j, a, y/) = -9(a, y/,j). (39) 

Thus the three 0 terms in (37) add to zero and we find that equation (25) also holds for a 
truncated Fourier representation. We can similarly multiply (35) by a (k) and (36) by 
o)*(k) and add the results, and separately multiply (36) by o*(k), and sum these over k to 
show that (27) and (29) are also valid in for truncated Fourier representations. Therefore, 
in the ideal case (v = T| = 0), E, He and A are also invariants in A>space, 
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£ = -L£* 2 [|y/(k)| 2 +|a(k)| 2 ] (40) 

in u 

a) H c =- L T Za(-k)cy(k) b) A = -^Z| «(k) | 2 . (41) 

IN u 2 N k 


(Note that each Zk sums independent coefficients twice.) The properties of the 0 
summations, as given by (38) and (39), ensure that E, He and A are the only quadratic 
invariants present in the truncated Fourier representation of 2-D homogeneous MHD 
turbulence. (If an ad hoc constant magnetic field were put into the 2-D problem, then A 
would no longer be conserved. We will not consider this situation further here.) 

5. MHD Turbulence as a dynamical system. In 2-D &-space, where k has integer 
components, there are a finite number of k such that |k| < A: max , and so there are a finite 
number of independent real and imaginary parts of the coefficients tf(k) and 0)(k). This 
number (call it Nr) is the dimension of the phase space T of the dynamical system 
associated with the truncated Fourier model of 2-D MHD turbulence. Since the equations 
of motion of the phase flow are (35) and (36), the divergence of the phase flow is 


a) 


1^2 = -V*2 

8a>(k) 


b) 


do(k) 

9a(k) 



(42) 


Clearly, for ideal MHD where v = r| = 0, the phase flow is divergenceless and the 
quantities E, He and A are conserved. The canonical distribution function D [19] for 
phase points satisfies a continuity equation 


dD dDaf k) dDdjk) 

dt u 3<o(k) 8o(k) 


= 0 . 


(43) 


Using equations (42) allows (43) to be written as 


dD 

dt 


: ir +z 

dt k 


®(k) 


dJD 

8cy(k) 


+ 8(k) 


dD 

8a(k) 



( 

\ 

= (v + ri) 

Zk 1 


V k 

/ 


D. 


(44) 


This tells us that if phase space is filled with points representing the initial states of 
systems [with values of E, He and A given by (40) and (41) depending on initial 
location], then the density of points in any volume moving with the phase flow increases 
exponentially with time, i.e., the gas of phase points collapses towards the origin, and the 
associated phase volume shrinks, when v and r| are positive definite constants. 
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6. Ideal MHD turbulence 

6.1 . Absolute equilibrium ensemble theory. When v = r| = 0, however, dD/dt = 0 and a 
Liouville theorem exists [21] Thus, since D is a constant function of the phase 
coordinates a(k) and 0)(k), it can only be a function of other phase functions which are 
themselves constant: D = D(E, He, A). Since E, He and A are additive, and since D must 
be multiplicative when the phase system T is split into parts T| and T 2 (i.e., D = D\Di) 
the form of D must be [19] 

a) D = Z'exp(-aE - p// c - y. A) b ) Z = jexp(-a E - (3 H c - y A)dT. (45) 

Here a, P and y are called inverse temperatures and the normalizing factor Z 1 is given in 
terms of the partition function Z. Equations (45) describe the statistical mechanics of a 
canonical ensemble [13, 19]. The partition function is a product of modal partition 
functions zr( k) and z/(k) and the volume is a product of modal volumes: 

a) Z = riZfl(k)z/(k) b) dT = n^T^OOdT, (k) (46) 

k k 

dT s (k) = da s (k )dco s (k), S = R,I (47) 

z/?(k) = z/(k) = Jexp[-oc(/f W +A 2 a 2 ) - Paw - ya 2 ] du>da. (48) 

The subscripts R, I in (48) denote real and imaginary parts, respectively. Expressions (40) 
and (41) for E, He and A are inserted in (45b) to produce (48), where the dummy 
variables of integration have been simplified. The integrations in (48) are performed 
between the limits — oo and +co for each a and w. The partition function thus depends on 
a, P, y, and k 2 . Carrying out the integrations, we find 

*g(k)= . ” S = R,I. (49) 

Va 2 -p 2 /4 + ay/* 2 

Let us define the expectation value of a phase quantity <|> as follows: 

<<])> = Z 1 J<|>exp(-a£- $H c -yA)dT. (50) 

Using the relations developed above, for 4) = w or a, it can be shown that 

a) <w(k)> = 0 b) <a(k)> = 0. (51) 

Canonical ensemble theory thus predicts that the mean values of the Fourier coefficients 
are zero. Second-order moments are also easily found to be 
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<lo>s( k)] 2 > = 2 S = R,I 

2(a 2 -p 2 /4 + ay/k 2 ) 

(52) 

([a s ( k)] 2 )= . ^ S = R,I 

2 (a 2 -p 2 IA + aylk 2 ) 

(53) 

(c°s (k) a s (k)> = — ~J~~ ry > S = R,I. 

a 2 -p 2 /4 + ay Ik 2 

(54) 


The values a, (3 and y are found by placing the expressions (52), (53) and (54) into (40) 
and (41), assuming E, He and A are constant and solving the resulting nonlinear 
equations. 

The ideal results presented above form what is called absolute equilibrium ensemble 
theory for 2-D MHD turbulence [10, 18]. (If an ad hoc uniform mean field B 0 is 
introduced into the 2-D case, A is no longer an invariant, although E and He remain so 
[30, 31].) In addition to these ideal 2-D MHD results, there are corresponding formulas 
for 3-D MHD turbulence [7, 32, 33, 37], as well as 2-D [17, 32, 33] and 3-D [16, 32, 33] 
Euler ( i.e ., ‘ideal Navier-Stokes’) turbulence. In 3-D MHD, the ideal invariants are E, He 
and the magnetic helicity [6, 23] Hm = 14[a-b] (if B 0 A 0, Hm is no longer invariant [33]). 
In 2-D homogeneous Euler turbulence, the integral invariants are E and £2 [17], while in 
3-D homogeneous Euler turbulence, E and the kinetic helicity [2, 23] Hk — !4[u-co] are 
invariant. The presence of at least one ideal invariant helicity in each case unites 2-D 
MHD with 3-D turbulence. 

Absolute equilibrium ensemble theory is a straightforward and useful application of 
Gibbsian statistical mechanics. The initial formulation of the absolute equilibrium theory, 
developed in the 1970s [7, 10, 16, 17] has been extended in several ways [30-35. These 
extensions and their uses are described in the following sections. 

6.2. Entropy. In the ideal case, if the expectation values of (40) and (41) are taken, and 
(52), (53) and (54) are used, three algebraic equations result: 

a <E> + $<H C > + Y <A> = -^r = n r (55) 

IN 2 

2a <H C > + (3<£ M > = 0 (56) 

a «£> - 2<E m » -y<A> = 0. (57) 

In the above, the expectation values <E>, <Hc> and <A> have essentially their initial 
values E, He and A because these have only minor canonical fluctuations. The value 
<Em> will, however fluctuate significantly since it is not an ideal invariant. Rename this 
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as (p- <Em>, for brevity. Then, solving the exactly determined equations (55), (56) and 
(57) for the unknown (p gives [30, 32, 33] 


a) a = 


n r <P 


<p(£-<p)-ff£ 


b) P = -2 


Hc_a 

<P 


c) y = — — —a. (58) 


Now, instead of having three variables to solve for, there is just one: cp. If the inverse 
temperatures in (58) are put into the distribution function D in (45), the result is D{tp), 
which is equal to the true distribution function when the equilibrium value of (p is found. 

The method for determining (p is to minimize the entropy functional o{tp) = - <ln D{<p)>, 
with respect to <p. Note that although entropy increases as previously isolated systems 
come together and interact, the entropy of an isolated system takes a minimum value with 
respect to temperature [13] (or temperatures, as in this case [34]). Equations (58) show 
that the inverse temperatures are really functions of only one variable, so that minimizing 
with respect to a, (3 and y is equivalent to minimizing with respect to (p. Using (50) and 
(58), along with (40) and (41) gives 

o{(p) = - <lnD(^?)> = « r (1 + A 2 ln7t)- 4-£ln(a 2 -{p 2 +ay k~ 2 ) (59) 

k 

Again, a, P and y are functions of <p as given in (58). Finding the value of (p for which 
minimizes o{<p) (there is only one minimum [36]) allows the equilibrium values of a, P 
and y to be found, and hence the expectation values (52), (53), and (54). This value is q> 
= <Em>, and the entropy of the dynamical system is S = g(<.Em>)- 

6.3. Broken ergodicity. The helicities He, Hm and Hk, mentioned above have one 
critical characteristic: they are pseudoscalars ( i.e ., they change sign) under either one or 
both of the symmetry transformations P {parity : x —> -x, so that u, a, j — > -u, -a, -j, but 
yr, b, 0) — > ys, b. to) and C {charge reversal : a, b, j -A -a, -b, -j). However, the equations 
of motion, such as (14) and (15) [or (35) and (36)], are unaffected by P or C. Also, 
although the distribution function D as given by (45) appears to be affected by P or C, it 
is not, for the inverse temperature P associated with He is clearly a pseudoscalar under P 
or C, as is seen in (58b), where P = - 2Hc<x/ (p (a and (p= Em are scalars, while He is a 
pseudoscalar). Similar relations can also be shown to exist for the invariant helicities in 
3-D ideal MHD and Euler turbulence [32, 33]. 

Thus, both the equations of motion and the statistical theory are invariant under P and C. 
However, in any solution of the ideal equations of motion the various invariant helicities 
always start with a definite sign (+ or -) and, because they are invariant, maintain that 
sign. In regard to ideal homogeneous turbulence as a dynamical system, each invariant 
helicity induces a disjoint structure on the available regions of phase space [35]. In ideal 
2-D MHD and 3-D MHD with B 0 * 0, as well as in 3-D Euler turbulence, there is one 
invariant helicity so that the available phase space is split into two effectively disjoint 
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parts. In 3-D MHD turbulence with B 0 = 0, there are two invariant helicities and the 
available phase space is effectively split into four disjoint parts. Only in the case of 2-D 
Euler turbulence, where there are no invariant helicities, is there a single, connected part 
in phase space for the system point to move on [32], 

The phrase effectively disjoint is used above because in a canonical ensemble, as we have 
here, invariant surfaces in phase space are actually ‘fuzzy,’ i.e., the phase point is not 
precisely confined the to invariant hypersurface T/ of dimension N, = Nr - 3 defined by 
the intersection of the Nr - 1 dimensional hypersurfaces of precisely constant E, He and 
A. (If phase motion were confined to T/ we would have a microcanonical ensemble [13, 
19].) Instead, E, He and A are canonically invariant , as they have very small fluctuations 
about a mean value, similar to the energy of a small physical system embedded in a large 
thermal heat bath (the ‘heat bath' here is the digital computer). The invariant 
hypersurface T/ is thus slightly spread out into an N r dimensional subspace AT in which 
the phase point resides with a probability given by the integral \^DdT = 1 , where D is the 
canonical distribution function given in (45). Since D is so highly peaked on T/, the 
probability that the phase point can actually be found between the effectively disjoint 
subspaces AT+ (AT = AT+ uAT-) associated with ±Hc is essentially zero. 

Numerical simulation has demonstrated the disjoint nature of phase space for ideal 
turbulence when invariant helicities are present [32, 33, 35]. Although (51) predicts that 
the mean values of co(k) and a(k) zero for 2-D MHD turbulence, in any dynamical 
simulation time averages of these modal values are no longer zero. The phase point 
begins and remains within one disjoint part of phase space, while ensemble averages are 
taken over all of phase space consistent with the canonically invariant values of E, \Hc\ 
and A. Thus, the symmetry of the governing equations and statistical theory under P and 
C is dynamically broken and the dynamical system is non-ergodic (as time averages do 
not match ensemble averages). This situation is often referred to as broken ergodicity 
[26]. Since the time-averaged values of co(k) and o(k) are no longer zero, these time 
averages define the existence of a coherent structure. 

The essential result is that direct numerical simulation of ideal, homogeneous turbulence 
is described by a dynamical system with a canonical distribution function. Furthermore, 
the presence of helical invariants induces broken ergodicity and coherent structure. Does 
this ideal turbulent phase space structre manifest itself in real turbulence when v and r| 
are non-zero? This is an important question that we will address presently. First, a brief 
discussion of real turbulence is given. 

7. Real MHD turbulence. In real 2-D MHD turbulence, when v ^ 0, r| 0, the 
spectrum is expected to be quite different from (52), (53) and (54). For forced, dissipative 
turbulence, a stationary state can be attained by injecting energy at some small value of k , 
after which it cascades through intermediate values of k with negligible dissipation (the 
inertial range), and is lost to heat at large k in the dissipation range, where k > ko, the 
dissipation wave number. (There can also be ‘inverse cascades’ to the smallest values of 
k.) The kinetic energy spectrum En{k) and magnetic energy spectrum E/u(k) are defined as 
integrals over azimuthal angle <() in &-space: 
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2n 


Efc(k) = -j J | u(k) | 2 kd<\> = — \ u)(k) \ 2 av 


(60) 


2 71 rrr 

EJ,k) = i 1 1 b(k) | 2 kd<S} = -\ y'(k) \ 2 ave 
o k 


(61) 


In the case of homogeneous Navier-Stokes turbulence, Kolmogorov [14] and Obukhov 
[24] used a dimensional arguments to predict that in the inertial range 

E K (k) = k- 513 . (62) 

For MHD turbulence with a mean field B 0 * 0, Iroshnikov [11] and Kraichnan [15] have 
predicted that Alfven wave effects produce inertial range spectra of the form 

E K (k) = Edk)~k~ V2 . (63) 

Here, we have no explicit mean field, since a large B 0 was invoked only to turn 3-D into 
2-D MHD turbulence, so that (62) should apply here, rather than (63). 

In the later stages of decay, theory [1] predicts a specific form for the energy spectrum 
Efc(k) of homogeneous Navier-Stokes turbulence. Since the equations for velocity and 
induction are the same when the non-linear terms are neglected, the magnetic energy 
spectrum E^k) of homogeneous MHD turbulence will take a similar form. The spectra 
for the later stages of decay (i.e., as / — > 0) are thus expected to behave as 

a) Efc(k) * A: 4 exp(-2v^ t) b) E^k) « k 4 exp(-2r\l^ t ). (64) 

Finally, the dissipation wave number ko for MHD turbulence can be approximated as 
(through dimensional reasoning similar to the Navier-Stokes case) 


k q — 271 


2 2 \ 

V T) ' 

V 2Q + T7 y 


- 1/4 


(65) 


The above results pertaining to both ideal and real turbulence can be used to examine 
some new ‘computer experiments.’ 

8. Numerical Simulation. We solve the £-space ordinary differential equations (35) and 
(36) for co(k) and a(k) with |k| < k max and reset all coefficients with |k| > k max to zero after 
each time step ( isotropic truncation). Let the number of k (with integer components not 
all zero) in the ball |k| < A: m ax be K. Associated with each of these k is one «(k) and one 
u>(k). Since a (k) = a(-k) and co (k) = co(-k), the number of wave vectors with 
independent coefficients is KJ 2. Since each a(k) or oo(k) has a real and an imaginary part. 
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the number of dimensions in the phase space is Nr = 4A72 = 2 K. Approximately, we have 
K = 7t^ m ax, although the exact number requires summation (this is Gauss’s circle 
problem), and numerical counting can arrive at precise values of K. 

Here, we present recent numerical simulations performed on an NxN grid of points with 
N = 512. Time integration is done by a third-order 'partially corrected Adams-Bashforth' 
scheme [8], and the linear dissipative terms are treated implicitly. Spatial differentiation 
is performed in &-space, where it is simple multiplication since the x-space V — > ik in k- 
space. Non-linear terms are evaluated by a) finding derivatives within cofactors in k- 
space, b) transforming these cofactors to x-space by fast Fourier transform (FFT), c) 
forming products at each point in x-space into Jacobians [such as J(a,i//) given by (16)] 
and d) transforming back to &-space by FFT. The aliasing that may occur in this 
transform method is fully eliminated in ideal runs by the Patterson-Orszag technique 
[27], where each evaluation of a non-linear term requires two sets of FFTs on shifted 
grids (this is a spectral method with ^max = 2 N 2 /9). The effects of aliasing are neglected 
in real runs (with only onset of FFTs per evaluation - this is a pseudospectral method 
with k 2 max = 'AN 2 -V2). Thus, for N= 512, ideal runs have k 2 max = 58254 and K= 183052 
(jikr m ax = 183010 - low by 42), while real runs have k 2 ^ = 65535 and K = 205856 
(nk 2 max = 205884 - high by 28). 

The simulation set consists of one ideal case and four real cases. All five runs had 
equivalent initial conditions (i.e., equal for k 2 < 58254). The dimension of the dynamical 
system in all cases is again Nr = 2 K. Thus, Nr = 366 1 04 in the ideal simulation, and Nr = 
41 1712 in the real simulations. Table I gives the values of v, r| and the magnetic Prandtl 
number Pm — v/r| for each of these five runs (the ideal case is denoted by Pm = ‘0/0’). 


Run 

V 


Pm 

11 

0 

0 

0/0 

R1 

0.001 

0.010 

0.10 

R2 

0.001 

0.004 

0.25 

R3 

0.0025 

0.0025 

1.00 

R4 

0.004 

0.001 

4.00 


Table I. Parameters for the numerical simulations 

The initial conditions were such that Ek = Em= 0.5, He = 0.024613, A = 0.031031 at / = 
0. The initial values of the complex coefficients co(k) and y(k) had random phase and 
magnitudes that varied as 

|co(k)| = |y(k)| ~ k 2 exp(-A^/16). (66) 

These initial values gave initial energy spectra Ex{k) - E/J^k) ~ k 6 exp(-A^/8), i.e., one 
that was highly peaked around k 2 = 5. The time step for numerical integration was A / = 
10‘ 4 , and each run in Table I was taken to / = 25. Each time step for the ideal run took 
about 9 cpu-seconds, while the time steps for the real runs took about 5 seconds each. 
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Thus, the total cpu-time invested (on a DEC Alpha running a 64-bit Fortran code named 
mhd2.f) was approximately 2000 cpu-hours and the run-time memory required was 66.4 
Mbytes. 

9. Numerical Results. First, the initial and time-averaged values, as well as standard 
deviations (as percentages of the time averages) of E, He, A and Em for the ideal run are 
given in Table 2. 


Integral 

Initial 

Time average 

Std. Dev. 

E 

1.0000 

1 .0004 

0.02 % 

He 

0.024613 

0.024603 

0.03 % 

A 

0.031031 

0.031031 

2x 1 O' 5 % 

Em 

0.50000 

0.56045 

6.35 % 


Table 2. Integral invariants for the ideal run (Pm = 0/0); averages from / = 0 to 25. 

Upon using the average values in Table 2 for E, He and A, a minimum for the entropy 
defined in equation (59) is found to occur at (p= <Em> = 0.515678. Using (58) then gives 
values to the inverse temperatures of a = 1.4371, P = -0.13713, and y = -1.4334. Putting 
these values into (52) and (53) allows us to predict the ideal spectra. In Figure 1, we 
compare these predicted spectra with spectra found in the ideal simulation, at times t = 1 
and t = 25. Simulation spectra were found by using (60) and (61), with discrete averaging 
done in bands around integer values of k n : n - Vi < k„ < n + 'A, where k„= 1,2, ..., [£ max ] 
(here, [ z ] is the integer part of the real number z). 

It is evident in Figure 1 that the spectra in the ideal run have evolved by / = 25 so as to be 
essentially the same as the predicted spectra at high k. However, they differ at low k, and 
the question is if and how quickly the low k spectra will collapse to the prediction. Let us 
look at Figure 2, where the ratio Ek/Em = E/Em- 1 versus time t is given for the five runs. 
For the ideal run (Pm = 0/0), this ratio does not appear to be approaching Em = 0.5 1 5678, 
i.e., Iogio(£x /Em) = -0.0272444, but rather is approaching logioCfix /Em) = -0.0492. 

The results shown in Figure 2 for Pm = 0/0 indicate that the low k coefficients may not 
evolve into their absolute equilibrium predictions. The value of Em appears to have 
become stationary at t = 20, and if we average from t = 20 to 25, we get the results in 
Table 3. As before, upon using the average values in Table 3 for E, He and A, a minimum 
for the entropy defined in equation (59) is found to occur at cp= <Em> = 0.515834. Using 
(58) then gives values to the inverse temperatures of a = 1.4366, p = -0.13699, and y = 
-1.4329. Putting these values into (52) and (53) allows us again to predict the ideal 
spectra, which does not perceptively differ from that shown in Figure 1. The predicted 
value <Em> ~ 0.515834 differs from the time average in Table 3 by more than 10 
standard deviations. Thus, there a small but significant departure in the simulation results 
from the predictions of the absolute equilibrium ensemble theory [7, 10, 16, 17]. 
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Figure 2 indicates that Ek/Em increase with decreasing P M , i.e., with increasing r|. Now, 
consider the total energy (26a) with respect to time, as shown in Figure 3. The ideal run 
(. Pm~ 0/0) has constant energy, as it should, while the real runs decay more quickly the 
lower their value of Pm- Thus decay is directly correlated with magnetic diffusivity r|, 
rather than with kinematic viscosity v, so that energy loss is primarily ohmic. (In the real 
runs, the dissipation wave number (59) was kp ~ k max at t ~ 0.5 and ko < &max otherwise.) 

The mismatch between absolute equilibrium ensemble theory and numerical simulation 
in ideal MHD turbulence is, again, caused by the broken ergodicity discussed in Section 
6.3. Broken ergodicity is clearly manifested in the time behavior of the coefficients a(k) 
and co(k) for ideal MHD turbulence. In order to see whether a similar effect appears in 
real MHD turbulence, we plot the evolution of o(k) and co(k) for several values of k. 

In Figure 4, the time evolution of o(k) for k = (0,1) is presented for all five runs. In 
Figure 5, the time evolution of ou(k), also for k = (0,1) is shown for all runs. These are the 
projections of the phase trajectory onto two-dimensional planes. It is clear that the phase 
trajectory is not even approximately symmetric about the origin. The effect is quite 
evident for k = 1 , though it diminishes for higher k. In Figure 6, the time evolution of a( k) 
(multiplied by k ) for k = (1,2) is presented for all five runs, where the projected phase 
trajectory does not seem as non-ergodic as for the k = 1 cases. In Figure 7, the time 
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evolution of co(k), also for k = (-2,2) is shown for all runs. In this last figure, the 
trajectory appears to flow into a more ergodic-looking shape around the origin. It must be 
remembered that these figures show different two-dimensional projections of a phase 
trajectory that evolves in a phase space of dimension of 366104 in the ideal run, and a 
dimension of 411712 in the real runs. The essential result is that the broken ergodicity 
that exists in ideal MHD turbulence also seems to appear in the real case. 



Integral 

Initial 

Time average 

Std. Dev. 

E 

1.0000 

1 .0007 

0.004 % 

He 

0.024613 

0.024593 

0.005 % 

A 

0.031031 

0.031031 

0.000 % 

E M 

0.50000 

0.53037 

0.252 % 


Table 3. Integral invariants for the ideal run (Pm = 0/0); averages from / = 20 to 25. 


Another result that appears in Figures 4 - 7 is that the low k modes grow in amplitude, 
while Figure 3 shows that total energy is dissipated in the real runs. Thus, we have an 
inverse cascade [10], in which energy flows to low k modes at the same time it is flowing 
to high k modes where it is lost. This can be quantified by plotting modal energy E(k) 
versus time, as is done in Figures 8 for E(k), k - 1; in Figure 9 for E^k), k = 1; and in 
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Figure 10, for E(k), k = 5. Figure 8 shows that the energy for k = 1 reaches a stationary 
state, while Figure 9 shows that the magnitude at k = 1 is mostly due to magnetic energy. 
Figure 10 illustrates the behavior of E(k) for k - 5 (and generally for E(k) with k > 1), 
which is to enter a stage of decay by t = 25. 



Figure 3. Total energy E with respect to simulation time. 


Figures 4-10 suggest that the coefficients r/(k) and co(k) are random variables with non- 
zero mean values. To examine this further, rolling time-averages and variances of the real 
and imaginary parts of both a(k) and co(k) were kept. Thus, for all coefficients, means 
and standard deviations were determined from t = 0 to t — 25. Non-zero mean values give 
rise to coherent structure, while variances quantify the random fluctuations within MHD 
turbulence. In the same way that complete coefficients determine kinetic and magnetic 
energy and spectra, the mean coefficients can be used to find coherent energy and 
spectra, and the variances can be used to find random energy and spectra. In Figure 1 1, 
coherent and random energy are plotted for the ideal case Pm = 0/0, and compared with 
absolute ensemble predictions. There is clear evidence that a significant amount of 
coherent energy exists across the whole spectral range, and particularly at low k. 
Although a more definitive result for this ideal case requires a considerable extension of 
the run with time, definitive results have been found previously in long runs on 32 2 grids, 
where significant amounts of coherent energy and structure were first discovered [32], 

The more interesting question at this point is the relevance of ideal results to real MHD 
turbulence. As Figures 4-10 show, there is a high degree of similarity at low k between 
the ideal and real runs presented here (where all runs began with essentially the same 
initial conditions). We now look at the final energy spectra of these five runs for further 
evidence of similar behavior. In Figure 12, we see the real and ideal total energy spectra 
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at t = 25, where it is clear that there is great similarity at low k and none at high k. Also, 
the Kolmogorov-Obukhov law [14,24] E(k ) ~ k' 5/i has been inserted for comparison. 



Figure 4. Evolution of tf(k), k = (0,1), for all runs. 



Figure 5. Evolution of co(k), k = (0,1), for all runs. 
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Figure 6. Evolution of a( k), k = (1,2), for all runs. 



Figure 7. Evolution of 0)(k), k = (-2,2), for all runs. 
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Figure 8. Evolution of E(k), k = 1 , for all runs. 



Time 

Figure 9. Evolution of E/J^k), k= 1, for all runs. 
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Time 


Figure 1 0. Evolution of E(k), k = 5, for all runs. 
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10. Conclusion. The representation of MHD turbulence by numerically realized finite 
dynamical systems produces intriguing results. The ‘broken ergodicity’ that exists in 
ideal MHD turbulence, and is strongest at low k, appears to have a counterpart in 
numerical simulations of real MHD turbulence. In fact, the low k behavior of the ideal 
and real runs presented here show a striking level of similarity. At high k, however, the 
energy spectra of ideal and real MHD turbulence differ radically. Thus, there is an 
apparent dichotomy in the in the transition from zero to non-zero values of v and r). First, 
ideal and real simulations agree well at low k (where ideal simulations do not agree with 
absolute ensemble theory due to broken ergodicity). Second, ideal and real simulations 
disagree at high k (where the ideal simulation agrees well with absolute ensemble theory). 
It seems as if the coefficients to(k) and n(k) at low k are only weakly coupled to the high 
k part of the spectrum, in both the ideal and real cases. This, in turn, provides a 
justification for the study of so-called 'large eddy simulations’ [29], in which the effects 
of high k coefficients are replaced by analytical sub-grid scale models. 

The result of varying the magnetic Prandtl number also leads to important insights. First, 
energy dissipation in MHD turbulence is primarily ohmic. Second, equipartion of energy 
between kinetic and magnetic modes in the real runs occurs at a value of Pm such that 
0.10 < Pm < 0.25, the precise value of which requires extended run times. Lastly, growth 
and saturation of the k = 1 inodes appears to depend only weakly on Pm- 

Although it is natural to extend simulations such as those presented here to significantly 
longer run times, the results found in the appear to be robust in light of the evident 
stationarity achieved during the latter part of these runs. There is much to be done and 
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there are also aspects of MHD turbulence other than those discussed here that can be 
profitably studied through numerical simulation. An additional benefit of 2-D MHD 
turbulence simulations is that these serve as representatives of 3-D MHD and Navier- 
Stokes turbulence, since these (but not 2-D Navier-Stokes turbulence) all possess at least 
one ideal invariant helicity, and that these numerical studies allow for higher & max than 
possible in 3-D simulations. 

While homogeneity is useful theoretical and computational assumptions, physical 
systems never truly support homogeneous flows. Nevertheless, The numerically observed 
growth and maintenance of coherent structures at the largest wavelengths is a qualitative 
feature which may represent a fundamental cause of the tendency of geophysical and 
astrophysical systems to spontaneously form large scale vortices or magnetic fields. 
Thus, broken ergodicity in MHD turbulence may well hold the key to understanding, for 
example, the ubiquitous presence of large scale magnetic fields in our earth, in our sun, 
and in other planets and stars. 
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